# broom package
#   tidy statistical models into data frames
library(broom)
library(ggplot2)
library(dplyr)

# 1. Linear regression ----

fit <- lm(mpg ~ wt + qsec, mtcars)
summary(fit)

# Statistics about each of the coeffieitns fit by th model
tidy(fit)
#   a new column 'term', means the data can be combined with other models

# per-observation info, such as fitted values and residuals
head(augment(fit))
#   add statistical info into original data

# compute per-model statisticals, such as R^2, AIC and BIC
glance(fit)

td <- tidy(fit, conf.int = TRUE)
ggplot(td, aes(estimate, term, color = term)) +
    geom_point() +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
    geom_vline(aes(xintercept = 0))

mtcars %>% group_by(am) %>% do(tidy(lm(mpg ~ wt, .)))

# 2. Visualization with ggplot2 ----
library(glmnet)
set.seed(03-19-2015)

# generated data with 5 real variables and 45 null, on 100 observations
nobs <- 100
nvar <- 50
real <- 5
x <- matrix(rnorm(nobs * nvar), nobs)
beta <- c(rnorm(real, 0, 1), rep(0, nvar - real))
y <- c(t(beta) %*% t(x)) + rnorm(nvar, sd = 3)
glmnet_fit <- cv.glmnet(x, y)

tidied_cv <- tidy(glmnet_fit)
glance_cv <- glance(glmnet_fit)

ggplot(tidied_cv, aes(lambda, estimate)) +
    geom_line(color = "red") +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
    scale_x_log10() +
    geom_vline(xintercept = glance_cv$lambda.min) +
    geom_vline(xintercept = glance_cv$lambda.1se, lty = 2)


library(survival)
surv_fit <- survfit(coxph(Surv(time, status) ~ age + sex, lung))
td <- tidy(surv_fit)
ggplot(td, aes(time, estimate)) +
    geom_line(size = 1) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2)
